Internal states of model isotropic granular packings. 
II. Compression and pressure cycles. 
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This is the second paper of a series of three investigating, by numerical means, the geometric 
and mechanical properties of spherical bead packings under isotropic stresses. We study the effects 
of varying the applied pressure P (from 1 or 10 kPa up to 100 MPa in the case of glass beads) 
on several types of configurations assembled by different procedures, as reported in the preceding 
paper 0]. As functions of P, we monitor changes in solid fraction "1>, coordination number z, 
proportion of rattlers (grains carrying no force) xo, the distribution of normal forces, the level of 
friction mobilization, and the distribution of near neighbor distances. Assuming that the contact 
law does not involve material plasticity or damage, <I> is found to vary very nearly reversibly with P 
in an isotropic compression cycle, but all other quantities, due to the frictional hysteresis of contact 
forces, change irreversibly. In particular, initial low P states with high coordination numbers lose 
many contacts in a compression cycle, and end up with values of z and xq close to those of the most 
poorly coordinated initial configurations. Proportional load variations which do not entail notable 
configuration changes can therefore nevertheless significantly affect contact networks of granular 
packings in quasistatic conditions. 

PACS numbers: 45.70.-n, 83.80.Fg, 46.65. +g, 62.20.Fe 
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I. INTRODUCTION 

The mechanical properties of solidlike granular pack- 
ings are traditionally studied, at the macroscopic level, 
in engineering fields such as soil mechanics |3j 
and are currently being investigated, with some atten- 
tion to the grain scale and micromcchanical origins of 
macroscopic behaviors, in condensed matter physics and 
material science communities [5, 6, 7]. 

The present paper, the second of a series of three, in- 
vestigates, by numerical simulations, the mechanical and 
microstructural response of a model material, the packing 
of identical spherical beads, to pressure intensity varia- 
tions. It refers a lot to the results of the previous, com- 
panion paper but may be read independently. 

Although molecular dynamics (or "discrete element") 
approaches have repeatedly been applied to sphere pack- 
ings d, [^, [13, UH in, EH, many important questions 
related to the microscopic origins of their macroscopic 
mechanical behavior in the quasistatic regime have not 
been fully explored yet. One such issue is the influence of 
the initial state, which is determined by the assembling 
process. In the first paper of the present series [2] (here- 
after referred to as paper I), the results of several pack- 
ing preparation methods, all producing ideally isotropic 
states, are compared. Direct compressions of granular 
gases produce states that do not depend on dynamical 
parameters if the compression is slow enough. Their solid 
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fraction $ and coordination number z* (evaluated on ex- 
cluding the rattlers, a proportion xq of grains which do 
not carry any force) are decreasing functions of the fric- 
tion coeflScient yU, from $ ~ 0.639 and z* = 6 for /i = 0, 
in which case the random close packing state (RCP) is 
obtained, down to $ ~ 0.593 and z* ~ 4.5 for ^ — 0.3. 
In paper I [l| we accurately checked the uniqueness of 
the RCP, on confronting our own numerical results with 
those of several recent publications, in which different 
numerical procedures were implemented [3, [isj . In the 
presence of intergranular friction, however, quite differ- 
ent packing states might be prepared. First, it is of 
course possible, in a simulation, to increase the friction 
coefficient once the packing is equilibrated under some 
pressure; such a numerical procedure can be regarded 
as a model for an assembling process in the presence 
of a lubricant within intergranular gaps in the labora- 
tory. Ideally, whatever the value of the friction coeffi- 
cient used to model the quasistatic mechanical proper- 
ties of the material, it is possible to assemble the sample 
with /i = (thus assuming ideal, perfect lubrication in 
the fabrication stage) and hence with the RCP density 
and coordination number. Once the grains are packed 
and form a solid material, contacts between grains can 
then be attributed the final, finite friction coefficient used 
in quasistatic modelling. Experimentally, it is of course 
well known that given granular materials can be packed 
with varying densities. A common method to make them 
denser, other then lubricating the contacts in the assem- 
bling stage, is the application of vibrations or "taps" . A 
numerical idealized vibration procedure, apt to prepare 
dense samples with little computation time, was defined 
in paper I. Surprisingly, although it produces isotropic 
states with densities close to the RCP value, their co- 
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ordination numbers are as low as in the loosest states 
assembled by direct compression. The small geometric 
differences between configurations with the same solid 
fraction but very different coordination numbers is still 
not accessible to tomographic observation techniques [l[ . 
Only mechanical properties can thus be confronted to 
experimental results, to determine whether or in which 
conditions the investigated numerical systems are close 
to experimental reality. 

Before studying elastic properties in paper III of the 
present series [1^1, one should first investigate the ef- 
fect of an isotropic compression. The application of a 
large enough confining pressure, usually at least a few 
tens of kPa (with rare exceptions [13, lll|)> is necessary 
before the macroscopic mechanical behavior of solidlike 
granular packings is tested 0, [l^, [l^], and characteris- 
tic quantities such as dilatancy and internal friction an- 
gle are measured. Experimental data on elastic mod- 
uli [m, m m, m, HB, 26, 27, I23 are also extremely 
scarce below that range. Most relevant laboratory sample 
histories to be understood in order to relate the macro- 
scopic response to internal variables and micromechan- 
ics involve an assembling stage, and then a compression 
stage, which is often isotropic or oedometric. It is there- 
fore necessary to assess the influence of pressure changes 
on the initial states. 

In addition, the material behavior under varying 
isotropic stress is interesting per se. The behavior of 
sands is traditionally regarded 0, |3, [2^ as elastoplastic 
under isotropic loading, with pressure cycles entailing ir- 
reversible density increases. Such effects are nevertheless 
considerably smaller than in cohesive materials such as 
clays @, H, |31j or powders [1^. It is worth investigat- 
ing such behavior in model sphere packings by numerical 
means. 



II. MODEL MATERIAL, MICROMECHANICAL 
PARAMETERS 



A. Contact model 

We briefly recall here the model material and the con- 
tact laws, which are described in paper I with more de- 
tails. Equal-sized spherical beads of diameter a (whose 
value, as we ignore gravity, will prove irrelevant), interact 
in their contacts by point forces of elastic, frictional and 
viscous origins. The Hertz law relates the normal elastic 
force N to the normal deflection h (approach of sphere 
centers closer than a) as : 



(1) 



E 



-, E being the Young mod- 



with the notation E = 

ulus of the beads, and ly the Poisson ratio. The Hertz 
law introduces a normal stiffness Kn — ^ that depends 
on h or on N. 



Tangential elasticity and friction are described with a 
simplified form of the Cattaneo-Mindlin-Deresiewicz re- 
sults [3^, in which the tangential stiffness Kt, relating 
the tangential elastic force increment to the relative tan- 
gential elastic displacement diiT in the contact, is pro- 
portional to K^: 

— (2) 



dT . , 2 

Kt — —. — = otKn with ax = 



The Coulomb condition with friction coefhcient fj, re- 
quires T to be projected back onto the circle of radius fiN 
in the tangential plane whenever the increment given by 
Eqn. ^ would cause its magnitude to exceed this limit. 
In order to avoid unphysical increases of elastic energy, 
T is scaled down in proportion with Kt when the elastic 
normal force N decreases, as indicated in paper I and 
advocated in [3l|. Tangential contact forces also move 
with the particles in contact, so that the condition of 
objectivity is satisfied (see paper I and ref. [s^). 

A viscous term opposing normal relative displacements 
reads (positive normal forces are conventionally repul- 
sive) : 

N'" = a{h)h, (3) 

with a damping coefficient a depending on elastic normal 
deflection h (or on elastic repulsive force TV), such that 
its value is a fixed fraction C of the critical damping co- 
efficient of the normal (linear) spring of stiffness K^ih) 
joining two beads of mass m: 



a{h) = Cy/2mKN{h). 



(4) 



We do not introduce any tangential viscous force, and im- 
pose the Coulomb inequality to elastic force components 
only. The main justification of such a term is computa- 
tional convenience (to accelerate the approach of equilib- 
rium states), and we could check that its value did not 
affect the statistical results on the configurations of the 
packings. 

The present numerical study was carried out with the 
elastic parameters E = 70GPa and v = 0.3 that are 
suitable for glass beads, and the friction coefficient is set 
to fi — 0.3. 



B. Stress control 

The numerical results presented below were obtained 
on samples of n = 4000 beads, enclosed in a cubic or par- 
allelipipedic cell with periodic boundary conditions. The 
sizes of the cell are denoted as La, parallel to coordinate 
axes a (1 < a < 3). Lq's vary simultaneously with the 
grain positions and orientations until mechanical equi- 
librium of all particles with the prescribed values Sq, of 
all three diagonal components aaa of the Cauchy stress 
tensor, 1 < a < 3, is obtained. One then has : 
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Here = L1L2L3 is the sample volume, 's are the 
coordinates of vector joining the center of bead i to 
the one of its contacting neighbor j (with the nearest im- 
age convention of periodic cells) and -P/^^'s are those of 
the corresponding contact force. This force is actually 
exerted by i onto j, so that the convention used is that 
tensile stresses are negative. Velocities of grain cen- 
ters comprise, in addition to a periodic field, an affine 
term corresponding to the global strain rate. Equations 
of motion for dimensions La are written in addition to 
the ordinary equations for the dynamics of a collection 
of solid objects, and they drive the system towards an 
equilibrium state in which condition ([5]) is obeyed. 

In the present study we always impose isotropic 
stresses, i.e. hydrostatic pressures P: Sq = P for a = 1, 
2, 3. 



C. Dimensionless parameters 

In addition to include friction coefficient /i and vis- 
cous dissipation parameter C, the important dimension- 
less control parameters for sphere packings under given 
pressure P are the reduced stiffness n and the inertia 
parameter /. k is chosen such that the typical contact 
deflection h is proportional to k^^. 
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(6) 



a correspondance which can be made accurate thanks to 
the relation 



P = 



(7) 



between pressure P = tra/S and the average normal 
force (N) in the contacts. ([7]) is exact provided h <^ a 
in all contacts and intercenter distances are taken equal 
to the diameter a. Here z denotes the cordination num- 
ber, equal to 2; = 2Nc/n, with the total number of 
force-carrying contacts in the packing. Rattlers, in pro- 
portion xq, have no such contact. We refer to te force- 
carrying network - the packing devoid of its rattlers - 
as the backbone, and to z*, which simply relates to z 
as rz = (1 — xo)z*, as the backbone coordination num- 
ber. Brackets denoting averages over all force-carrying 
contacts, one has 
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The limit of rigid grains is approached as k ^ 00. 

K can be used to determine whether the material within 
the grains is likely to be imposed stresses beyond its elas- 
tic limit. The maximum pressure, at the center of a 
Hertzian contact between spheres of diameter a, carrying 
a normal force iV, is [301 



2 X 31/3 ^2/3 
TT a2/3 



iVl/3. 



Under pressure P, corresponding to k by ([6]), when the 
average normal force in contacts is (N), one can deduce 
from ([7|) 
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7r2/3{z$)i/3 l^(7V)y ■ ^ ' 



Likewise, the maximum shear stress Tmax, which is 
reached inside the grains near the contact region will 
be m (for 1/ = 0.3) 



E 



E 



(9) 



Eqns. [5]and[n]show that very high stress levels, up to a 
non-negligible fraction of elastic modulus E are reached 
if K, is not large enough. With our choice of material 
parameters for glass beads, we get ~ 0.051 for 

P=10 MPa and k-^/^ ~ 0.11 for P=100 MPa, while the 
numerical prefactor is only slightly lower than 1 ('^ 0.8) 
if z$ = 4 (a typical value) in ([5]). Such high stresses are 
very likely to entail particle breakage or plastic strains 
(according to the materials the grains are made of). 

In our simulations we set our lowest pressure level for 
the simulation of glass beads to 1 kPa or 10 kPa, corre- 
sponding to K ~ 181000 and k ~ 39000 with the elastic 
properties of glass. This enables us to explore the entire 
experimental pressure range, and to approach the large k 
limit too. Up to the maximum pressure value 100 MPa, 
we assume elastic contact behavior, but one should be 
careful on comparing the numerical results in the higher 
pressure states [P > 10 MPa) to experimental ones. 

Dynamical effects are assessed on comparing the strain 
rate e to intrinsic inertial times, such as the time needed 
for a particle of mass m, initially at rest, accelerated by 
a typical force Pa^, to move on a distance a. This leads 
to the definition of a dimensionless inertia parameter : 



I — e\J m/aP. 



(10) 



The quasistatic limit can be defined as / 0. / is a con- 
venient parameter to describe internal states and write 
down constitutive laws f or g ranular materials in dense 
shear flow [H [H, HE [H, [33| . 



D. Initial states 

The present paper is devoted to the study of the influ- 
ence of quasistatic pressure changes to granular packings 
assembled by different means, as described in paper I [l|. 
Four different states were prepared under low pressure, 
and some of their basic characteristics are recalled in ta- 
ble m Such state variables are monitored in the follow- 
ing as a function of pressure in isotropic compression or 
pressure cycles. In addition to solid fraction propor- 
tion of rattlers xq, backbone (or force-carying structure) 
coordination number z* , Table U provides some global in- 
formation on force distributions. Z{2) is characteristic of 
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TABLE I: Isotropic states (k ~ 39000 for A and C, k ~ 181000 for B and D) for different assembling procedures. 



Procedure 


$ z- xo (%) Z(2) Ah Mi 


A 


0.6370 ± 0.0002 6.074 ± 0.0015 1.3 ± 0.2 1.53 


B (^io = 0.02) 


0.6271 ± 0.0002 5.80 ±0.007 1.95 ± 0.02 1.52 0.016 0.018 


C (vibration) 


0.635 ±0.002 4.56 ±0.03 13.3 ±0.5 1.65 0.135 0.181 


D 


0.5923 ± 0.0006 4.546 ± 0.009 11.1 ±0.4 1.58 0.160 0.217 



the width of the distribution of normal forces: 



Z(2) 



{If) 



(11) 



Ml and M2 are the average levels of friction mobilization 



{i.e., 



N 



-) for contacts carrying normal forces, respec- 



tively, larger and smaller than the average {N) . 

In paper I we also recorded other geometric data, in 
particular pair correlation functions and distributions of 
near neighbor gaps h. The latter can be expressed as gap- 
dependent coordination numbers, defining z{h) as the av- 
erage number of neighboring beads around a central one, 
separated by an interstice smaller than h. z{0) thus co- 
incides with the contact coordination number. Due to 
the rattlers, the proportion of which -see table |l]- can 
exceed 10% of the total number of grains, such geomet- 
ric data are however somewhat ambiguously defined: the 
positions of the rattlers are not fixed by the rigid back- 
bone. Thus one may define z^{h), on using the arbitrary 
positions obtained at the end of the simulation, when the 
packing first equilibrates within the prescribed numeri- 
cal tolerance. One then has 2:^(0) ~ z (recall z counts 
only force-carrying contacts) if the equilibrium state is 
accurately computed, because there are very few con- 
tacts bearing a normal force below tolerance. In an at- 
tempt to define more intrinsic geometric data, we defined 
z^^{h) in paper I 0] as the gap-dependent coordination 
number in the configuration obtained once all rattlers are 
pushed against the backbone, in random directions. In 
their new position, the rattlers now have three contacts 
with the backbone (except in the rare case when inter- 
rattler contacts are obtained). It was argued in paper 
I that the resulting structure was likely to resemble, to 
some extent, granular assemblies under gravity, when the 
weight of the grains is very small in comparison to the 
local stress. z^^(O) can be regarded as a geometric defi- 
nition of a contact coordination number (it is, in general, 
slightly larger than z* ^ z/{l — xq)). 



III. NUMERICAL RESULTS 

We first specify the numerical compression proce- 
dure in paragraph IIII Al then describe the effects of an 
isotropic compression and a pressure cycle in terms of 
global variables (Section IIII Bp as well as local geome- 
try (Section IIII C[) . We then test the simplest prediction 
scheme for the evolution of coordination number, that 



of homogeneous strain at the microscopic level, in Sec- 
tion ElDl 



A. Numerical procedure 

The results presented below pertain to equilibrium 
configurations at variable isotropic pressure P, obtained 
by a stepwise compression (respectively: decompression) 
process in which P, within the controlled stress scheme 
described in Section III Bl is increased (respectively: de- 
creased) by a factor vTO. In each pressure step a condi- 
tion of slow enough strain rate was enforced, so that the 
inertia parameter, as defined by (|10p with the currently 
imposed pressure level, was kept below a maximum value: 
/ < 10"'^ for compression, / < 10"'' for decompression. 
Such values were chosen to ensure independence of the 
results on dynamical parameters / and C- It was observed 
that a decompression process requested more care, due 
to its greater instability. Whereas a compression of the 
sample beyond its equilibrium density will be strongly 
opposed, at growing P, by elastic forces in the network, 
too large an expansion, as P decreases, might cause the 
contact network to break apart, resulting in a dynam- 
ical process similar to assembling a granular gas, when 
the externally applied pressure finally drives the system 
back to a denser equilibrium configuration. Such events 
might entail a significant remoulding of the contact net- 
work and large departures from equilibrium conditions. 
This should of course be avoided in a procedure designed 
to model a quasistatic evolution, as close as possible to 
the limit of small strain rates. 

Configurations are deemed equilibrated when, defin- 
ing £p = 10~^Pa^ as a small tolerance on forces and 
es = 10~''Pa'^ as a small tolerance on energies, the four 
following conditions are simultaneously satisfied : 

• each coordinate of the total force on each grain is 
smaller than ep', 

• each coordinate of the total moment on each grain 
is smaller than e^a; 

• all stresses have their prescribed values with a rel- 
ative error smaller than ep: 



{a = 1, 2, 3) 



PI 



P 



the kinetic energy per grain is smaller than ep- 
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To distinguish between the backbone and the rattlers, 
the same method is appHed as presented in paper I [l| . 

Such procedures were apphed to samples A to D below, 
with P ranging from its smallest value 1 kPa (for B and 
D, corresponding too k ~ 181000), or 10 kPa (for A and 
C, corresponding to k ~ 39000), up to 100 MPa {n ~ 84), 
and then back to its initial low value. Letters A, B, C, 
D will hereafter denote pressure-dependent configuration 
series. Although initial states A and B were assembled 
with coefficients of friction lower than the chosen value 
fi = 0.3, we study quasistatic compressions with fi = 0.3 
for all sample series. We regard the smaller friction levels 
applied to configurations A and B in the assembling stage 
as models for lubricated grains, and assume that the lu- 
bricant ceases to operate once solid particles finally touch 
one another, as in equilibrated packings and during qua- 
sistatic compression tests. As a reference for comparisons 
with other states, and because it was studied in the lit- 
erature [13, , we also prepared another configuration 
series we denote as AO, obtained from the initial A state 
on compressing a frictionless system (thus series A and 
AO share the same initial low-pressure state, but differ as 
soon as P is altered). 

All results are averaged over 5 samples of n=4000 
beads, and error bars correspond to one standard de- 
viation. 



B. Evolution of global state variables 

Figs. [T] and [2] display the evolution of solid fraction $, 
backbone coordination number z*, and rattler fraction 
xq in sample series A, B C and D in the pressure cycle. 
Fig. [H shows that the solid fraction change with pressure 
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FIG. 1: (Color online) Evolution of packing fraction as a func- 
tion of pressure P in glass bead packings (bottom axis), or 
dimensionless stiffness parameter (top axis), in (from top 
to bottom) states A (red crosses, continuous line), C (black 
square dots, continuous line) B (blue asterisks, dotted line) 
and D (green open squares, dotted line). 



is almost perfectly reversible: the data points correspond- 
ing to the compression and decompression parts of the 
pressure cycle are almost indistinguishable. More pre- 
cisely, once the pressure had returned to its lowest value 
in samples A to C, the packing fraction was observed to 
have changed by very small amounts, below 2 • 10^**. The 
loosest state, D, undergoes a slight compaction. Yet, this 
effect apparently decreased as the maximum prescribed 
value for parameter / was changed from 10"** to 10~^ 
upon unloading (the reported results corresponding to 
this latter value). Our model material thus differs from 
sands, which are reported to respond to such cycles with 
notable irreversible density increases 0, [1] ■ It should be 
noted, though, that we are using a contact model without 
plasticity or particle damage, which, as argued on eval- 
uating, in Sec. Ill C[ the maximum pressure and shear 
stress in the grains near contact points with Eqns. ([8]) 
and is quite unrealistic for the highest pressure lev- 
els simulated. Stress concentrations in contacts between 
angular particles like sand grains, with corners or asper- 
ities [lO, , are more severe than between smooth ob- 
jects and should enhance the effects of anelastic material 
behavior within the grains. The smallness of irreversible 
compaction in our simulations suggests that such macro- 
scopic behavior, in sands, originates in contact mechanics 
rather than in collective effects. 

The reversibility of the response to the pressure cycle is 
however only apparent, as the coordination number does 
not return to its initial value. 

As expected, z* increases under a growing confining 



pressure (Fig. 2(a) I : as the particles pack more closely in 
a smaller volume, near neighbors come into contact, z* 
reaches about 7.3 at the highest pressure in the densest 
samples, A and C. Correlatively, an increasing number 
of rattlers get trapped as their free volume shrinks, and 
are recruited by the force-carrying network. The initially 
large fraction of rattlers in states C and D {xq > 10%) 
steadily decreases as P grows( Fig. |2(b)[ ) and has virtu- 
ally disappeared at P = 100 MPa. 

The evolution of coordination numbers on unloading 
is more surprising. While low coordination states C and 
D exhibit a very limited hysteresis effect and eventually 
retrieve their initial, low z* values (about 4.6), with a 
slightly lower rattler fraction, samples of types A and B, 
in which z* was initially high, lose contacts as a result 
of the pressure cycle and end up with z* values below 5 
(about 4.8 for A, and 4.5 for B), closer to C and D ones 
than to where they started, with a substantial rise in the 
population of rattlers. (Let us recall that samples A and 
B are regarded in the study of quasistatic compression as 
made of frictional beads with /i = 0.3, like the others). 
The behavior of (frictionless) samples AO is of course dif- 
ferent, for they cannot be stable at low pressure below 
z* = 6 40]. Fig. [3] compares the evolutions of z* in A 
and AO series, and shows that z* is very nearly reversible 
in the AO series. The unloading curves in A states start- 
ing at lower pressures, 3.16 Mpa and 1 MPa instead of 
100 MPa, also shown on Fig. [31 witness a lower, but sig- 
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FIG. 2: (Color online) Backbone coordination number z* (a) 
and proportion of rattlers xq (b) as functions of P or , 
same symbols as on Fig. [1] 



nificant decrease of z* from its initial value z* ~ 6 at the 
end of the cycle. The shape of the force distribution and 
the mobilization of friction also change with P, as shown 
by the evolution of parameters Z{2), Mi, M2 on Fig. ID 
As a general rule, the width of the force distribution cor- 
relates with the level of force indeterminacy, relatively 
to the number of degrees of freedom. Contact elasticity 
tends to share forces rather evenly, because contact force 
values should minimize the intergranular elastic energy, 
subject to the constraint that they balance the applied 
pressure (this elastic energy as a function of forces is 
written further below in connection with a discussion of 
irreversibility in pressure cycles, and the minimization 
property is exploited in paper III [l^ to estimate bulk 
moduli). More precisely, the increments of forces due 
to pressure increases will tend to reduce the width of the 
distribution, the faster the less constrained the minimiza- 
tion, i.e. the larger the degree of force indeterminacy. 
Thus in configurations A, the large coordination number 




P (kPa) 

FIG. 3: z* versus P or in pressure cycle in series A 
(crosses) and AO (dots), showing reversibility for AO. Shorter 
cycles (up to 0.316 MPa and 1 MPa) than the one of Fig. [5] 
are also shown for A. 



enables a quick narrowing of the distribution under grow- 
ing pressure. In states C, the same tendency is present, 
but the evolution is much slower, as there are less pos- 
sibilities to distribute forces in a more tenuous network 
while maintaining equilibrium. However, C samples gain 
contacts faster than D ones (Fig. 2(a) I, for which the nar- 
rowing effect is even slower. Finally, the extreme case is 
the situation of isostaticity, as in the AO series, in which 
the distribution of forces is geometrically determined in 
the rigid limit of k — s- -f 00. As, furthermore, the increase 
of z with P is not very fast in that case, since z is already 
large from the beginning, the shape of the distribution 
remains nearly constant. A few normal force probabil- 
ity distribution functions at different pressure levels are 
shown on Fig. [5l 

The evolution of force values and friction mobilization 
on unloading is more complicated: all three parameters 
shown on Fig. [4] first increase, then go through a max- 
imum and end up, at the initial pressure value, with a 
value comparable with the initial one (except for fric- 
tion mobilization parameters Mi and M2 in A systems, 
because they started at zero). In a granular sample con- 
trolled in displacements or strains, rather than stresses, 
large self balanced forces can in some situations remain 
when the external load that created them is removed, the 
simplest ex amp le being that of one particle wedged in a 
corner [4l|, l42j |. Our observations indicate that such a 
phenomenon does not take place in a situation of con- 
trolled stress state: all forces are of the order of the 
average force, which is related to the current pressure 
by ([7]) , even though contacts have carried forces that were 
larger by orders of magnitude in the past. This suggests 
that the set of admissible contact forces, restricted to 
the intersection of an /i-dimensional affine space (due to 
equilibrium relations) with a cone (due to Coulomb in- 
equalities) is bounded. Yet during unloading many more 
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FIG. 4: (Color online) From bottom to top; Z{2), Ml and 
M2 versus P or in compression cycle. Symbols as on 
Fig. [T] for states A, C, and D. Series AO represented with 
(red) dots joined by dotted line for Z{2). Hysteresis loops for 
Z{2) first decrease, then increase back on unloading and go 
through a maximum (except for AO, in which cas it is nearly 
constant). Mi and M2 behave in a similar way, with the 
special circumstance that their initial values are equal to zero 
in A states (assembled without friction). 



sliding contacts are observed than at growing pressure, 
due to the effects of decreasing normal force components, 
and the level of friction mobilization is higher (Fig. [4]). 
Meanwhile, the distribution of normal forces gets wider. 
The global influence of the past loading, with contacts 
previously carrying larger forces, enhances force hetero- 
geneities. A related quantity is the elastic energy stored 
in the contacts. The total elastic energy per grain w 
reads (from Eqns. ([T|) and ^) 



1 " 

EE 



=1 



Once adimensionalized by Ea^, we denote it as w. On 
exploiting Eqn. ([7|) it is conveniently expressed as: 



32/3^5/3 z{5/3) 



(12) 



5 ^2/3(j)5/3^5/2 ■ 

In (fT^ . Z(5/3), related to force moments, is close to 



FIG. 5: From bottom to top, evolution of normalized force 
distributions P{f), with / = N/{N), with growing pressure 
in samples A, C, D, AO. P value in kPa are 10 (except for D: 
P = 1), 100, 10^, 10* and 10^. AU four distributions tend to 
narrow as P grows, but at very different rates. 



Z(5/3), which can be defined on replacing exponent 2 by 
5/3 in Eqn. (fTTj) . with the following slight modification. 
With ax defined in ^ as the constant ratio of tangential 
to normal stiffnesses, and with the notation ttn for the 



ratio 



N 



in a contact, let us define 



Z(5/3) = 



(7V5/3(1 



5ri 



■)) 



(iV)5/3 



(13) 



Z(5/3) thus depends on the force distribution and also 
on friction mobilization, although for fi — 0.3 its rela- 
tive difference with Z(5/3) is small (of the order of -1^, 
with Ml as plotted on Fig. [5]). The energy per parti- 
cle, w, scales as which is expected since this is 
proportional to h^^"^ for h (x the typical normal con- 
tact deflection, w is larger for low coordination num- 
bers (weaker networks) , and larger force disorder (higher 
Z(5/3)). (It should be recalled that we use pressure, 
rather than strain, as the control parameter, hence a 
larger elastic energy for softer materials). Thus in A 
configurations, w is larger, for given k, on decompressing, 
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P (fca) 



FIG. 6: (Color online) Increment of packing fraction A^"^"^ 
gained between the two states of equal pressure, reached at 
growing and at decreasing P, in states A (red, crosses) and C 
(black, square dots). Note the scale of density changes (A$ 
of order 10""). 



another manifestation of the irreversibility of the cycle. 
If we assume that the curve P($) is quasistatically fol- 
lowed up to the maximum pressure, and then exactly re- 
traced back on decompressing, this leads to a paradox, as 
some elastic energy appears to be gained at no expense. 
Thus one has to account for very small irreversible den- 
sity changes, for energetic consistency. Such changes in 
$, between the growing and decreasing pressure parts of 
the cycle are shown on Fig. [HI In the case of A con- 
figurations, one even observes a slight decompaction on 
decreasing P back to its lowest, initial value. Although 
surprising, this phenomenon should be expected in the 
rigid limit P — s- or k — s- oo, because as explained in 
paper I, the initial A configuration, which was assembled 
without friction, is a local maximum of $ subject to im- 
penetrability constraints. Another conclusion of paper 
I [l| is that the only way to increase density in such a 
sample is to produce, by enduring agitation or repeated 
shakes, notable traces of crystalline order. This should 
not happen in a slow, quasistatic compression experiment 
with only one pressure cycle. To check for energetic con- 
sistency, one may note on Fig. [5] however that the change 
of (f> is positive at high pressure. The total energy fed 
into the system in the cycle is 



Aw''''' 



TT r A$"' (P)dP 

6 J 1$^ ' 



(14) 



the integral running over the whole pressure interval of 
the compression cycle. Consequently (see Fig.|6]) the con- 
tribution of the irreversible increase of is largely dom- 
inant, because it is integrated over a much wider pres- 
sure interval. The small changes in density between the 
compression and the decompression curves at the same 
pressure values are large enough to explain the change in 
elastic energy, and that of potential energy as well when 



the cycle ends up decreasing the density (which happens 
for A samples). 



C. Pair correlations and near neighbor distances 

The smallness or absence of irreversible compaction in 
the pressure cycle implies that the samples do not avoid 
contact deflections by finding denser packing arrange- 
ments. Thus interparticle correlation patterns should 
witness favored near neighbor distances which typically 
scale like <I'^^/'^. This is shown for C configurations 
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(a)g(r) versus r/a in C configurations at difli'erent P. 
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FIG. 7: Pair correlation functions at P=10, 100, 1000, 
10'*, 10^ kPa in configurations C at growing pressure, with- 
out (top), and with (bottom) rescaling distance r as r* = 
r($/$o)'''^ 

on Fig. [71 on rescaling the distance axis, using coor- 
dinate r* — r($/$o)^^^ with $o the initial low pres- 
sure solid fraction, the different g(r) curves are super- 
imposed. In agreement with the observations made in 
paper I [l[, where the relationships between pair corre- 
lation functions and contact networks were discussed, a 
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closer look on such correlations will reveal differences in 
the details of the peaks associated with changes in the 
coordination number with <I>. Figs. [8] and [9] respectively 
show functions (h) and z^^ih) at growing P values, 
using the corresponding change of scale for interstice ft,, 
h* = ($/$o)^/^(a + h) - a. Those data suggest that 



show that such predictions, albeit reasonable, are not ex- 
act. In particular, the gradual capture of rattlers by the 




0.01 , . . 0.02 

h /a 



0.03 



FIG. 8: (Color online) Gap-dependent neigbor coordination 
number versus rescaled interstice h* = ($/4'o)^''^(a + 

ft) — a at different P (same as on Fig. O in states A (red) , C 
(black) and D ( green). 




0.03 



FIG. 9: (Color online) Same as Fig. [8]for definition z''{h) of 
the gap-dependent neigbor coordination number. 

the homogeneous shrinking of distances implied by the 
rescaling of abscissae on the graphs of Figs. |7(b)[ [8] and 
[9] is an approximation with some discrepancies at small 
intergranular distances. Curves corresponding to pres- 
sures other than the lowest one on Figs. [5] and O start at 
distance [($/$o)^/^ — l]a > and the corresponding val- 
ues of z(h) on the curve for the lowest pressure value are 
the predictions for the coordination number on assum- 
ing homogeneous shrinking strains. Differences therefore 



force-carrying network as P grows (see Fig. 2(b) ) cannot 



be adequately described by the homogeneous shrinking 
assumption: the rattlers will not start carrying forces 
when one interstice with a backbone grain is closed. The 
use of definition z^^ {h) should in principle improve this 
kind of prediction: once positioned against the backbone 
(with 3 contacts), the rattlers are much more likely to 
create new contacts bearing nonzero forces when they 
touch new neighbors. Yet, the improvement of curve su- 
perpositions on Fig. [5] compared to Fig. [5] is marginal. 
This suggests that the inaccuracy of the prediction of 
coordination numbers is not only due to the capture of 
rattlers by the growing backbone, but also stems from 
the failure of the assumption of homogeneous shrinking. 



D. Can one predict the changes in coordination 
number ? 



The results of the prediction of the coordination num- 
ber, assuming all distances uniformy shrink, are shown 
on Fig. [in] for systems A and C under growing pressure. 
The agreement is very good in state A (except at high 
pressure, where z is slightly underestimated), and fair in 
state C. For C configurations, the prediction was done 
separately for both z and z-^^(O), showing a somewhat 
better accuracy at low pressure in the second case. Unfor- 
tunately, the mechanically important coordination num- 
ber is 2^(0) = z rather than z^^(O). To evaluate P as a 
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FIG. 10: (Color online) Predictions for z — z^(0) in samples A 
and C, and for z^' (fi) in samples C, based on the homogeneous 
shrinking assumption. 

function of $, one needs to account for two phenomena: 
the increase of the elastic normal deflection in the con- 
tacts that already existed at the lowest pressure, and the 
creation of new contacts due to the closing of open inter- 
stices. Both effects are evaluated with the assumption of 
homogeneous rescaling of all distances according to the 
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density change, respectively exploiting the previous mea- 
surements of the distribution of sphere overlaps (related 
to that of normal forces), and of the function z{h) (with 
no significant difference in accuracy on using or z^^). 
The predicted values of z, although not very accurate for 
small changes of at low pressures, globally capture the 
marked growing trend above 1 MPa. The predictions of 
density increases are compared with the simulation re- 
sults on Fig. [nl showing good agreement (with a slight 
underestimation at high pressure). The prediction of P 
is understandably more accurate than that of the coordi- 
nation number, because it is not very sensitive, at first, 
to errors in the estimation of the density of newly created 
contacts, which initially carry very small forces. 




.62 0.63 0.64 0.65 0.66 0.67 



0.68 
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FIG. 11: (Color online) $ versus P or in samples A (red) 
and C (black). Dots: measurements. Dotted lines: predic- 
tions, based on the homogeneous shrinking assumption from 
the initial state of lowest pressure. 



FIG. 12: (Color online) Coordination number z versus "1> at 
decreasing P in samples A (red) and C (black). Dots: mea- 
surements. Dotted lines: predictions, based on the homoge- 
neous expansion assumption from the initial state of highest 
pressure. 
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One may also attempt to predict the decrease of coor- 
dination number in the decompression part of the pres- 
sure cycle. Such a prediction is based on the distribu- 
tion of particle overlaps (or contact deflections), rather 
than near neighbor distances. The relevant information 
is therefore the normal force histogram for the highest 
pressure level, as shown, e.g. on Fig. O However, this 
is a rather crude approximation, which leads to large er- 
rors for the coordination number variation with density, 
as shown on Fig. I12| and very poor predictions indeed 
for the coordination number relationship to the decreas- 
ing pressure, as apparent on Fig. 1131 Such an assumption 
of homogeneous expansion proves in particular unable to 
provide a correct estimate of the properties at low density 
or pressure, as it ignores the requirement of mechanical 
rigidity. We are not aware of a simple prediction scheme 
that would be able to provide a reasonably accurate de- 
scription the reduction of coordination number in the A 
state on reducing the confining pressure. 



FIG. 13: (Color online) Coordination number z versus de- 
creasing P (or in samples A (red) and C (black). Dots: 
measurements. Dotted lines: predictions, based on the homo- 
geneous expansion assumption from the initial state of highest 
pressure. 

IV. DISCUSSION 

The effect of a compression on the four series of 
isotropic packings we have been studying can be broadly 
summarized as the closing of additional contacts and the 
gradual reduction of the characteristic disorder of gran- 
ular systems, as witnessed by the narrowing of the force 
distribution (Figs. S] and [5]). Geometric changes con- 
form to the homogeneous shrinking assumption on large 
scale, and the resulting predictions for the near-neighbor 
distances and the coordination numbers are reasonable, 
if not very accurate, approximations (Figs. [TUl and [TT|l . 
even though they cannot correctly account for the re- 
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cruitment of rattlers (Fig. |2(b)[ ) by the growing back- 
bone. It proves difficult to accurately estimate small z* 
increases, to which, as will be studied in [l^ (paper III), 
shear moduli of poorly coordinated packings are espe- 
cially sensitive. The changes in the forces and the mo- 
bilization of friction are not appropriately described by 
such a simple model. On assessing the performance of 
the homogeneous shrinking approximation, one thus re- 
trieves the classification of length scales introduced in 
paper I P, Section IV.E.2]. Global chan ges on scales 
above about 0.05a appear to abide by the homogeneous 
strain assumption, hence the superposition of pair corre- 
lation functions on Fig. |7(b)[ Pair correlations between 
neighbors at smaller distances (or details of the peaks 
of g{r)) are only approximately predicted on rescaling 
all distances by the same factor (as appears on Figs. [5] 
and[9]). And small distances of the order of (contact 
deflections related to forces) do not abide by this homo- 
geneity of strain. Otherwise, on rescaling coordinates by 
a factor 1 — e, where <C e ^ 1, one would replace any 
contact deflection hhy ta + h, which for e ^ would 
result in a much stronger narrowing of the force distribu- 
tion than the one observed. This assumption of homo- 
geneous strain (or afEne displacements) will be further 
tested on dealing with elastic moduli in paper III [31 • 

The effects of a pressure reduction are more surpris- 
ing. Although the evolution of solid fractions departs 
very little from reversibility (Figs. [T]and[6]), large ini- 
tial coordination numbers in configurations A and B do 



not survive a pressure cycle (see Figs. 2(a) and [3]). Such 



effects are not predicted by the simple assumption of ho- 
mogeneous expansion, which grossly fails to reproduce 
the evolution of coordination number and density on re- 
ducing the confining pressure (Figs. [T^ and [T31). The 
memory of larger stresses, upon decompressing, imparts 
wider force distributions and larger friction mobilizations 
in some pressure range (Fig. 3]) , while such reductions of 
coordination numbers take place. It should be expected 
that decompression is less predictible, because it is an 
evolution towards a larger disorder, and small differences 
can be amplified in the process. This contrasts with the 
compression phase, in which, for instance, the differences 



between configurations A and C tend to disappear. Den- 
sity differences are recovered on decreasing P, with the 
additional phenomenon that new internal states at low 
pressure are thus being prepared, which also differ from 
the initially assembled ones. While this phenomenon es- 
capes the currently available modelling schemes, it can be 
noted that configurations with a high coordination num- 
ber, for nearly rigid grains (low pressure or high stiffness 
parameter k), are extremely rare, since each contact re- 
quires a new equation to be satisfied by the set of sphere 
centre positions. Equilibrium states of rigid, frictionless 
sphere assemblies, which are the initial states for config- 
uration series A, apart from the motion of the scarce rat- 
tlers, are isolated points in configuration space, because 
of isostaticity, as discussed in paper I [l|. As the pres- 
sure cycle, at the microscopic scale, is not reversible, due 
to friction and to geometric changes, one should not ex- 
pect such exceptional configurations to be retrieved upon 
decreasing the pressure. 

We thus conclude that the internal state of granular 
packings, in addition to the assembling process, the ef- 
fect of which was studied in paper I [l|, varies accord- 
ing to the history of stress intensities, even though, un- 
like in cohesive materials [1^ , and in contrast with 
changes in stress directions, such loading modes only en- 
tail very small irreversible strains. Such commonly used 
characteristics of granular packings as coordination num- 
ber, force distribution and friction mobilization level are 
sensitively affected by their compression history, while 
strains and density changes remain very small after the 
assembling stage. In particular, large coordination num- 
bers associated with an ideally successful suppression of 
friction in the sample preparation stage seem even more 
unlikely to occur generally in isotropic sphere assemblies 
close to the RCP density, because they do not survive 
compression cycles. Elastic properties are studied in pa- 
per III [l6j . where we relate them to the microstucture 
of such states, thereby allowing for compararisons of nu- 
merical results to experimental ones. 

As possible developments of the present study, one may 
simulate the effects of irreversible contact deformation, 
due to material plasticity or particle breakage. 
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